//============================================================================//
//                                                                            //
// TetGen                                                                     //
//                                                                            //
// A Quality Tetrahedral Mesh Generator and A 3D Delaunay Triangulator        //
//                                                                            //
// Version 1.6.0                                                              //
// August 31, 2020                                                            //
//                                                                            //
// Copyright (C) 2002--2020                                                   //
//                                                                            //
// Hang Si                                                                    //
// Research Group: Numerical Mathematics and Scientific Computing             //
// Weierstrass Institute for Applied Analysis and Stochastics (WIAS)          //
// Mohrenstr. 39, 10117 Berlin, Germany                                       //
// si@wias-berlin.de                                                          //
//                                                                            //
// TetGen is a tetrahedral mesh generator. It creates 3d triangulations of    //
//   polyhedral domains. It generates meshes with well-shaped elements whose  //
//   sizes are adapted to the geometric features or user-provided sizing      //
//   functions. It has applications in various applications in scientific     //
//   computing, such as computer graphics (CG), computer-aided design (CAD),  //
//   geometry processing (parametrizations and computer animation), and       //
//   physical simulations (finite element analysis).                          //
//                                                                            //
// TetGen computes (weighted) Delaunay triangulations for three-dimensional   //
//   (weighted) point sets, and constrained Delaunay triangulations for       //
//   three-dimensional polyhedral domains. In the latter case, input edges    //
//   and triangles can be completely preserved in the output meshes.  TetGen  //
//   can refine or coarsen an existing mesh to result in good quality and     //
//   size-adapted mesh according to the geometric features and user-defined   //
//   mesh sizing functions.                                                   //
//                                                                            //
// TetGen implements theoretically proven algorithms for computing the        //
//   Delaunay and constrained Delaunay tetrahedralizations. TetGen achieves   //
//   robustness and efficiency by using advanced techniques in computational  //
//   geometry. A technical paper describes the algorithms and methods         //
//   implemented in TetGen is available in ACM-TOMS, Hang Si ``TetGen, a      //
//   Delaunay-Based Quality Tetrahedral Mesh Generator", ACM Transactions on  //
//   Mathematical Software, February 2015, https://doi.org/10.1145/2629697.   //
//                                                                            //
// TetGen is freely available through the website: http://www.tetgen.org.     //
//   It may be copied, modified, and redistributed for non-commercial use.    //
//   Please consult the file LICENSE for the detailed copyright notices.      //
//                                                                            //
//============================================================================//

#ifndef tetgen_IO_H
#define tetgen_IO_H

// TetGen default uses the double-precision (64 bit) for a real number.
//   Alternatively, one can use the single-precision (32 bit) 'float' if the
//   memory is limited.

#define REAL double  // #define REAL float

// The maximum number of characters in a file name (including the null).

#define FILENAMESIZE 1024

// The maximum number of chars in a line read from a file (including the null).

#define INPUTLINESIZE 2048

// C standard libraries to perform Input/output operations, general utililities,
//   manipulate strings and arrays, compute common mathematical operations,
//   get date and time information.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <string>

// The types 'intptr_t' and 'uintptr_t' are signed and unsigned integer types,
//   respectively. They are guaranteed to be the same width as a pointer.
//   They are defined in <stdint.h> by the C99 Standard.

#include <stdint.h>
#include "tetgenbehavior.h"

// fastlookup() -- A fast, unsafe operation. Return the pointer to the object
//   with a given index.  Note: The object's block must have been allocated,
//   i.e., by the function newindex().
#define fastlookup(pool, index) \
  (void *) ((pool)->toparray[(index) >> (pool)->log2objectsperblock] + \
            ((index) & (pool)->objectsperblockmark) * (pool)->objectbytes)

class tetgenio
{
public:
    // A "polygon" describes a simple polygon (no holes). It is not necessarily
    //   convex. Each polygon contains a number of corners (points) and the same
    //   number of sides (edges).  The points of the polygon must be given in
    //   either counterclockwise or clockwise order and they form a ring, so 
    //   every two consecutive points forms an edge of the polygon.
    struct polygon {
        int* pVertexList{ nullptr };
        int NumVertex{ 0 };
    };

    // A "facet" describes a polygonal region possibly with holes, edges, and 
    //   points floating in it.  Each facet consists of a list of polygons and
    //   a list of hole points (which lie strictly inside holes).
    struct facet {
        polygon* pPolygonList{ nullptr };
        int NumPolygon{ 0 };
        REAL *pHoleList{ nullptr };
        int NumHole{ 0 };
    };

    // A "voroedge" is an edge of the Voronoi diagram. It corresponds to a
    //   Delaunay face.  Each voroedge is either a line segment connecting
    //   two Voronoi vertices or a ray starting from a Voronoi vertex to an
    //   "infinite vertex".  'v1' and 'v2' are two indices pointing to the
    //   list of Voronoi vertices. 'v1' must be non-negative, while 'v2' may
    //   be -1 if it is a ray, in this case, the unit normal of this ray is
    //   given in 'vnormal'. 
    struct voroedge {
        int v1, v2;
        REAL vnormal[3];
    };

    // A "vorofacet" is an facet of the Voronoi diagram. It corresponds to a
    //   Delaunay edge.  Each Voronoi facet is a convex polygon formed by a
    //   list of Voronoi edges, it may not be closed.  'c1' and 'c2' are two
    //   indices pointing into the list of Voronoi cells, i.e., the two cells
    //   share this facet.  'elist' is an array of indices pointing into the
    //   list of Voronoi edges, 'elist[0]' saves the number of Voronoi edges
    //   (including rays) of this facet.
    struct vorofacet {
        int c1, c2;
        int *elist;
    };

    // Additional parameters associated with an input (or mesh) vertex.
    //   These informations are provided by CAD libraries. 
    struct pointparam {
        REAL uv[2];
        int tag;
        int type; // 0, 1, or 2.
    };

    // Callback functions for meshing PSCs.
    typedef REAL (* GetVertexParamOnEdge)(void*, int, int);
    typedef void (* GetSteinerOnEdge)(void*, int, REAL, REAL*);
    typedef void (* GetVertexParamOnFace)(void*, int, int, REAL*);
    typedef void (* GetEdgeSteinerParamOnFace)(void*, int, REAL, int, REAL*);
    typedef void (* GetSteinerOnFace)(void*, int, REAL*, REAL*);

    // A callback function for mesh refinement.
    typedef bool (*TetSizeFunc)(REAL*, REAL*, REAL*, REAL*, REAL*, REAL);

    // Items are numbered starting from 'FirstNumber' (0 or 1), default is 0.
    int FirstNumber{ 0 };

    // Dimension of the mesh (2 or 3), default is 3.
    int iMeshDim{ 3 };

    // Does the lines in .node file contain index or not, default is 1.
    int iUseIndx{ 1 };

    // 'pPointList':  An array of point coordinates.  The first point's x
    //   coordinate is at index [0] and its y coordinate at index [1], its
    //   z coordinate is at index [2], followed by the coordinates of the
    //   remaining points.  Each point occupies three REALs. 
    // 'pPointAttributeList':  An array of point attributes.  Each point's
    //   attributes occupy 'NumPointAttribute' REALs.
    // 'pPointMtrList': An array of metric tensors at points. Each point's
    //   tensor occupies 'numberofpointmtr' REALs.
    // 'pPointMarkerList':  An array of point markers; one integer per point.
    // 'point2tetlist': An array of tetrahedra indices; one integer per point.
    REAL *pPointList{ nullptr };
    REAL *pPointAttributeList{ nullptr };
    REAL *pPointMtrList{ nullptr };
    int  *pPointMarkerList{ nullptr };
    int  *point2tetlist{ nullptr };
    pointparam *pPointParamList{ nullptr };
    int NumPoint{ 0 };
    int NumPointAttribute{ 0 };
    int numberofpointmtrs{ 0 };
 
    // 'pTetrList':  An array of tetrahedron corners.  The first 
    //   tetrahedron's first corner is at index [0], followed by its other 
    //   corners, followed by six nodes on the edges of the tetrahedron if the
    //   second order option (-o2) is applied. Each tetrahedron occupies
    //   'NumCorner' ints.  The second order nodes are ouput only. 
    // 'pTetrAttributeList':  An array of tetrahedron attributes.  Each
    //   tetrahedron's attributes occupy 'NumTetrAttribute' REALs.
    // 'pTetrVolumeList':  An array of constraints, i.e. tetrahedron's
    //   volume; one REAL per element.  Input only.
    // 'neighborlist':  An array of tetrahedron neighbors; 4 ints per element. 
    // 'tet2facelist':  An array of tetrahedron face indices; 4 ints per element.
    // 'tet2edgelist':  An array of tetrahedron edge indices; 6 ints per element.
    int  *pTetrList{ nullptr };
    REAL *pTetrAttributeList{ nullptr };
    REAL *pTetrVolumeList{ nullptr };
    int  *neighborlist{ nullptr };
    int  *tet2facelist{ nullptr };
    int  *tet2edgelist{ nullptr };
    int NumTetr{ 0 };
    int NumCorner{ 4 };
    int NumTetrAttribute{ 1 };

    // 'pFacetList':  An array of facets.  Each entry is a structure of facet.
    // 'pFacetMarkerList':  An array of facet markers; one int per facet.
    facet *pFacetList;
    int *pFacetMarkerList;
    int NumFacet{ 0 };

    // 'pHoleList':  An array of holes (in volume).  Each hole is given by a
    //   seed (point) which lies strictly inside it. The first seed's x, y and z
    //   coordinates are at indices [0], [1] and [2], followed by the
    //   remaining seeds.  Three REALs per hole. 
    REAL *pHoleList;
    int NumHole{ 0 };

    // 'pRegionList': An array of regions (subdomains).  Each region is given by
    //   a seed (point) which lies strictly inside it. The first seed's x, y and
    //   z coordinates are at indices [0], [1] and [2], followed by the regional
    //   attribute at index [3], followed by the maximum volume at index [4]. 
    //   Five REALs per region.
    // Note that each regional attribute is used only if you select the 'A'
    //   switch, and each volume constraint is used only if you select the
    //   'a' switch (with no number following).
    REAL *pRegionList;
    int NumRegion{ 0 };

    // 'refine_elem_list': An array of tetrahedra to be refined.  The first 
    //   tetrahedron's first corner is at index [0], followed by its other 
    //   corners. Four integers per element.
    // 'refine_elem_vol_list':  An array of constraints, i.e. tetrahedron's
    //   volume; one REAL per element.
    int  *refine_elem_list;
    REAL *refine_elem_vol_list;
    int  NumRefineElem{ 0 };

    // 'pFacetConstraintList':  An array of facet constraints.  Each constraint
    //   specifies a maximum area bound on the subfaces of that facet.  The
    //   first facet constraint is given by a facet marker at index [0] and its
    //   maximum area bound at index [1], followed by the remaining facet con-
    //   straints. Two REALs per facet constraint.  Note: the facet marker is
    //   actually an integer.
    REAL *pFacetConstraintList;
    int NumFacetConstraint{ 0 };

    // 'pSegmentConstraintList': An array of segment constraints. Each constraint 
    //   specifies a maximum length bound on the subsegments of that segment.
    //   The first constraint is given by the two endpoints of the segment at
    //   index [0] and [1], and the maximum length bound at index [2], followed
    //   by the remaining segment constraints.  Three REALs per constraint. 
    //   Note the segment endpoints are actually integers.
    REAL *pSegmentConstraintList;
    int NumSegmentConstraint{ 0 };

    // 'trifacelist':  An array of face (triangle) corners.  The first face's
    //   three corners are at indices [0], [1] and [2], followed by the remaining
    //   faces.  Three ints per face.
    // 'trifacemarkerlist':  An array of face markers; one int per face.
    // 'o2facelist':  An array of second order nodes (on the edges) of the face.
    //   It is output only if the second order option (-o2) is applied. The
    //   first face's three second order nodes are at [0], [1], and [2],
    //   followed by the remaining faces.  Three ints per face.
    // 'face2tetlist':  An array of tetrahedra indices; 2 ints per face.
    // 'face2edgelist':  An array of edge indices; 3 ints per face.
    int *trifacelist;
    int *trifacemarkerlist;
    int *o2facelist;
    int *face2tetlist;
    int *face2edgelist;
    int NumTriaFace{ 0 };

    // 'edgelist':  An array of edge endpoints.  The first edge's endpoints
    //   are at indices [0] and [1], followed by the remaining edges.
    //   Two ints per edge.
    // 'edgemarkerlist':  An array of edge markers; one int per edge.
    // 'o2edgelist':  An array of midpoints of edges. It is output only if the
    //   second order option (-o2) is applied. One int per edge.
    // 'edge2tetlist':  An array of tetrahedra indices.  One int per edge.
    int *edgelist;
    int *edgemarkerlist;
    int *o2edgelist;
    int *edge2tetlist;
    int NumEdge{ 0 };

    // 'vpointlist':  An array of Voronoi vertex coordinates (like pPointList).
    // 'vedgelist':  An array of Voronoi edges.  Each entry is a 'voroedge'.
    // 'vfacetlist':  An array of Voronoi facets. Each entry is a 'vorofacet'.
    // 'vcelllist':  An array of Voronoi cells.  Each entry is an array of
    //   indices pointing into 'vfacetlist'. The 0th entry is used to store
    //   the length of this array.
    REAL *vpointlist;
    voroedge *vedgelist;
    vorofacet *vfacetlist;
    int **vcelllist;
    int numberofvpoints{ 0 };
    int numberofvedges{ 0 };
    int numberofvfacets{ 0 };
    int numberofvcells{ 0 };

    // Variable (and callback functions) for meshing PSCs.
    void *geomhandle;
    GetVertexParamOnEdge getvertexparamonedge;
    GetSteinerOnEdge getsteineronedge;
    GetVertexParamOnFace getvertexparamonface;
    GetEdgeSteinerParamOnFace getedgesteinerparamonface;
    GetSteinerOnFace getsteineronface;

    // A callback function.
    TetSizeFunc tetunsuitable;

    // Input & output routines.
    bool LoadTetrMesh(std::string &sName);
    bool load_node_call(FILE* infile, int markers, int uvflag, char*);
    bool load_node(char*);
    bool load_edge(char*);
    bool load_face(char*);
    bool load_tet(char*);
    bool load_vol(char*);
    bool load_var(char*);
    bool load_mtr(char*);
    bool load_elem(char*);
    bool load_poly(char*);
    bool load_ply(char*);
    bool load_plc(char*, int);
    bool load_tetmesh(char*, int);
    void save_nodes(const char*);
    void save_elements(const char*);
    void save_faces(const char*);
    void save_edges(char*);
    void save_neighbors(char*);
    void save_poly(const char*);
    void save_faces2smesh(char*);

    // Read line and parse string functions.
    char *readline(char* string, FILE* infile, int *linenumber);
    char *findnextfield(char* string);
    char *readnumberline(char* string, FILE* infile, char* infilename);
    char *findnextnumber(char* string);

    static void init(polygon* p) {
        p->pVertexList = (int *) NULL;
        p->NumVertex = 0;
    }

    static void init(facet* f) {
        f->pPolygonList = (polygon *) NULL;
        f->NumPolygon = 0;
        f->pHoleList = (REAL *) NULL;
        f->NumHole = 0;
    }

    // Initialize routine.
    void initialize()
    {
        pPointList = (REAL *) NULL;
        pPointAttributeList = (REAL *) NULL;
        pPointMtrList = (REAL *) NULL;
        pPointMarkerList = (int *) NULL;
	    point2tetlist = (int *) NULL;
        pPointParamList = (pointparam *) NULL;
    
        pTetrList = (int *) NULL;
        pTetrAttributeList = (REAL *) NULL;
        pTetrVolumeList = (REAL *) NULL;
        neighborlist = (int *) NULL;
	    tet2facelist = (int *) NULL;
	    tet2edgelist = (int *) NULL;
    
        trifacelist = (int *) NULL;
        trifacemarkerlist = (int *) NULL;
        o2facelist = (int *) NULL;
        face2tetlist = (int *) NULL;
	    face2edgelist = (int *) NULL;
    
        edgelist = (int *) NULL;
        edgemarkerlist = (int *) NULL;
        o2edgelist = (int *) NULL;
        edge2tetlist = (int *) NULL;
    
        pFacetList = (facet *) NULL;
        pFacetMarkerList = (int *) NULL;
    
        pHoleList = (REAL *) NULL;
        pRegionList = (REAL *) NULL;
    
        refine_elem_list = (int *) NULL;
        refine_elem_vol_list = (REAL *) NULL;
    
        pFacetConstraintList = (REAL *) NULL;
        pSegmentConstraintList = (REAL *) NULL;
    
        vpointlist = (REAL *) NULL;
        vedgelist = (voroedge *) NULL;
        vfacetlist = (vorofacet *) NULL;
        vcelllist = (int **) NULL;

        tetunsuitable = NULL;

        geomhandle = NULL;
        getvertexparamonedge = NULL;
        getsteineronedge = NULL;
        getvertexparamonface = NULL;
        getedgesteinerparamonface = NULL;
        getsteineronface = NULL;
    }

    // Free the memory allocated in 'tetgenio'.  Note that it assumes that the 
    //   memory was allocated by the "new" operator (C++).
    void clean_memory()
  {
    int i, j;

    if (pPointList != (REAL *) NULL) {
      delete [] pPointList;
    }
    if (pPointAttributeList != (REAL *) NULL) {
      delete [] pPointAttributeList;
    }
    if (pPointMtrList != (REAL *) NULL) {
      delete [] pPointMtrList;
    }
    if (pPointMarkerList != (int *) NULL) {
      delete [] pPointMarkerList;
    }
	if (point2tetlist != (int *) NULL) {
      delete [] point2tetlist;
    }
    if (pPointParamList != (pointparam *) NULL) {
      delete [] pPointParamList;
    }

    if (pTetrList != (int *) NULL) {
      delete [] pTetrList;
    }
    if (pTetrAttributeList != (REAL *) NULL) {
      delete [] pTetrAttributeList;
    }
    if (pTetrVolumeList != (REAL *) NULL) {
      delete [] pTetrVolumeList;
    }
    if (neighborlist != (int *) NULL) {
      delete [] neighborlist;
    }
    if (tet2facelist != (int *) NULL) {
	  delete [] tet2facelist;
	}
	if (tet2edgelist != (int *) NULL) {
	  delete [] tet2edgelist;
	}

    if (trifacelist != (int *) NULL) {
      delete [] trifacelist;
    }
    if (trifacemarkerlist != (int *) NULL) {
      delete [] trifacemarkerlist;
    }
    if (o2facelist != (int *) NULL) {
      delete [] o2facelist;
    }
    if (face2tetlist != (int *) NULL) {
      delete [] face2tetlist;
    }
	if (face2edgelist != (int *) NULL) {
      delete [] face2edgelist;
    }

    if (edgelist != (int *) NULL) {
      delete [] edgelist;
    }
    if (edgemarkerlist != (int *) NULL) {
      delete [] edgemarkerlist;
    }
    if (o2edgelist != (int *) NULL) {
      delete [] o2edgelist;
    }
    if (edge2tetlist != (int *) NULL) {
      delete [] edge2tetlist;
    }

    if (pFacetList != (facet *) NULL) {
      facet *f;
      polygon *p;
      for (i = 0; i < NumFacet; i++) {
        f = &pFacetList[i];
        for (j = 0; j < f->NumPolygon; j++) {
          p = &f->pPolygonList[j];
          delete [] p->pVertexList;
        }
        delete [] f->pPolygonList;
        if (f->pHoleList != (REAL *) NULL) {
          delete [] f->pHoleList;
        }
      }
      delete [] pFacetList;
    }
    if (pFacetMarkerList != (int *) NULL) {
      delete [] pFacetMarkerList;
    }

    if (pHoleList != (REAL *) NULL) {
      delete [] pHoleList;
    }
    if (pRegionList != (REAL *) NULL) {
      delete [] pRegionList;
    }
    if (refine_elem_list != (int *) NULL) {
      delete [] refine_elem_list;
      if (refine_elem_vol_list != (REAL *) NULL) {
        delete [] refine_elem_vol_list;
      }
    }
    if (pFacetConstraintList != (REAL *) NULL) {
      delete [] pFacetConstraintList;
    }
    if (pSegmentConstraintList != (REAL *) NULL) {
      delete [] pSegmentConstraintList;
    }
    if (vpointlist != (REAL *) NULL) {
      delete [] vpointlist;
    }
    if (vedgelist != (voroedge *) NULL) {
      delete [] vedgelist;
    }
    if (vfacetlist != (vorofacet *) NULL) {
      for (i = 0; i < numberofvfacets; i++) {
        delete [] vfacetlist[i].elist;
      }
      delete [] vfacetlist;
    }
    if (vcelllist != (int **) NULL) {
      for (i = 0; i < numberofvcells; i++) {
        delete [] vcelllist[i];
      }
      delete [] vcelllist;
    }
  }

    // Constructor & destructor.
    tetgenio() {initialize();}
    ~tetgenio() {clean_memory();}

}; // class tetgenio

#endif // #ifndef tetgenH
